	subroutine cd_init
      include 'cst.inc'
c*********************************************************************************
c	  this subroutine is only called once at the beginning of the
c	  model run. its purpose is to provide a depth-, and time-
c	  averaged entrained sediment concentration (ccstar) for the
c	  whole range of depths in the model grid.  this is then held
c	  constant in all subsequent time-steps because subroutine
c	  candd is used to compute cstar, dl and dc as the sea floor
c	  evolves.
c	     
c	  sets the large-scale diffusion coefficients (dlbar) &
c       (dcbar) and the time-averaged areal entrained sediment
c       concentration (ccstar) for all of the grid. in the initial
c       program dlbar=f(y),dcbar=g(y) and ccstar=fo(h).
c       l and c are the local longshore and cross-shore
c       directions determined from the average longshore
c       current direction. the ccstar function was determined
c       in earlier numerical experiments to yield an offshore
c       depth profile like the one used in program gridset.
c
c       Output variables of this subroutine: ccstar,hcs,jlimit,dc,dl
c********************************************************************************	
	real cs(imax,jmax)
c
c     Ddetermine the cross-, and long-shore large-scale diffusion coefs.
c     The calculation below is based equation (4) of Niedoroda(1995), although 
c     the equation is for 1-D case. It assumes that the diffusion coefficient
c     is a function a distance offshore or along shore 
c
	do i=1,imax                       !Loop along the shore (x-direction)
	   do j=1,jshore(i)               !Loop from the domain boundary to the shoreline node
	     yy=(jshore(i)-j)*deltay    !distance from the shoreline
	     dcc=dczero+dcyyy*yy        !dc is for cross-shore, dcyyy is defined in definits and a read-in parameter 
	     dc(i,j)=min(dcmax,dcc)     !truncate the dc value to a specified maximum
	     dll=dlzero+dlyyy*yy        !dl is for long-shore
           dl(i,j)=min(dlmax,dll)     !truncate the dl value to a specified maximum
        enddo
      enddo
c
c    determine cstar = f(h) using function that yields an equilibrium
c    bypass shelf profile of the form used in setgrid.  using this with
c    a grid developed from setgrid will make the initial profile close 
c    to a equilibrium bypassing shape.
c
c	Beta - Scales the exponent in the wave-drift
c	Wo - Scales the wave-drift sed. trans.
c	DC(I,J) - Cross-shore diff. coef.in flow coords.{1} 
c	H(I,J,iTime) - Depths in grid, fill index in surf-zone cells {SG}
c	Cstar(I,J,iT) - time-,& depth-averaged entrained sed.concentration {1}
c
c     The calculation is based on Equation (18) of Niedoroda (1995)
c     Since this is for 2_D case, dcdy here is equivalent to dcdx in Niedoroda (1995)
c
  	do i=1,imax
	   nn=jshore(i)-1
	   nj=jshore(i)
	   cs(i,nj)=0.	
   	   do j=nn,1,-1
	      dcdy=-wo*(exp(beta*((h(i,j)+h(i,j+1))/2.))/  !h here is negative 
     +                         ((dc(i,j)+dc(i,j+1))/2.))
	      cs(i,j)=dcdy*deltay+cs(i,j+1)
	   enddo
	enddo      
c	Hmax - Max.(ie. most negative) depth in the surf zone cell (SG)
c	Hmin - Min. depth in the surf zone cell	 (SG)
c??????????????? The hmax is different from that of SetGrid
c
c --- Find the maximum depth of the entire grid ---
c
	hmax=-1.
	do i=1,imax
	   nnn=jshore(i)
	   do j=1,nnn
	      if(h(i,j).lt.hmax) then
	        hmax=h(i,j)
	        ihmax=i           !The role # of hmax
	        jlimit=jshore(i)  !The column # of hmax 
	      endif
	   enddo
	enddo        
c
c --- Assign the depth of the role of hmax to hcs ---
c
	nnn=jshore(ihmax)
      do j=1,nnn
	   hcs(j)=h(ihmax,j)
	   ccstar(j)=cs(ihmax,j)
	enddo   
c	
c --- add one more really deep point to file (to be sure that all depths are covered)
c?????????????? Why shifting two points here
	nnn=nnn+1
	do j=nnn,2,-1
	   hcs(j)=hcs(j-1)
	   ccstar(j) = ccstar(j-1)
	enddo   	
	hcs(1)=2*hcs(2) !double the depth
	ccstar(1)=ccstar(2)
c
c--- add one more point near zero water level to be sure all points are covered
c
	fac = (ccstar(nnn)-ccstar(nnn-1))/(hcs(nnn)-hcs(nnn-1))
	nnn=nnn+1
 	hcs(nnn)=3.0
	ccstar(nnn)=ccstar(nnn-1)+fac*(hcs(nnn)-hcs(nnn-1))

	jlimit=nnn

	cmin = 99999.
	do j=1,nnn
      	 cmin=amin1(ccstar(j),cmin)
      enddo	 
	do j=1,nnn
	   ccstar(j)=ccstar(j)-cmin
      enddo

c	do 72 j=1,nnn
c72	write(*,*)j,hcs(j),ccstar(j)
c	pause

	return
	end